home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / multi.i < prev    next >
Text File  |  1995-12-12  |  49KB  |  1,309 lines
  1. /*
  2.    MULTI.I
  3.    Drat/DSP routines defining and implementing a standard for
  4.    retrieving opacity/emissivity data from multiple files and
  5.    combining it for a single problem.
  6.  
  7.    $Id$
  8.  */
  9.  
  10. /*
  11.    Design for a Drat interface which handles opacity/emissivity
  12.    input from several post-processing files.  Key features:
  13.  
  14.    (1) One file is the "master file", which contains all of the
  15.        zones (rt, zt, ireg from the original calculation), at all
  16.        of the times that were dumped.  This file may or may not
  17.        contain its own gb, ekap, and akap arrays.
  18.  
  19.    (2) All other files may contain only subsets of the zones or
  20.        times of the master file.  If they contain only a subset of
  21.        the times, the transport calculation will be done only at
  22.        the times which are common to every file.  Any zone may
  23.        get contributions to its opacity and emissivity from zero,
  24.        one, or more of the files.
  25.  
  26.    (3) Every file (except possibly the master file) contains its own
  27.        gb, ekap, and akap arrays.  Optionally, it may contain a gexist
  28.        array, which is one shorter than gb, and is 0 in any bins which
  29.        do not exist in that file.  The group dimension of ekap and akap
  30.        should be the number of non-zero elements of gexist.  This allows
  31.        calculations with non-contiguous coverage of photon energy.
  32.        Alternatively, codes which have no concept of photon energy bins
  33.        may write a gav array, and omit the gb array.  A suitable gb will
  34.        be inferred, and gexist is meaningless for this case.
  35.        ===>> In no case may gb or gexist (or gav) change either as a
  36.              function of time, or from zone to zone.
  37.        However, the group structures of the various files need not match.
  38.        You can either specify an overall group structure, or have the
  39.        code generate one for you before you perform the transport
  40.        calculation; in either case, the code will interpolate from
  41.        the group structure in each file onto the common group structure
  42.        where the transport calculation is performed.
  43.  
  44.    (4) Zones are determined by a separate list of "master" zone indices
  45.        for each slave file (and optionally for the master file, too).
  46.        If a zone in some slave file is marked in its zone list, that
  47.        file will contribute both opacity and emissivity to the total
  48.        for that zone, otherwise it will not contribute.
  49.  
  50.    (5) The quantities currently dumped, akap (1/length) and
  51.        ekap (power/photon energy/solid angle/area) may not be appropriate
  52.        for another opacity code, because the ekaps from several sources
  53.        cannot simply be added.  Therefore, the slave files may optionally
  54.        dump the variable emiss (power/photon energy/solid angle/volume)
  55.        instead of ekap.  The emiss and akap from several different files
  56.        simply add.
  57.  
  58.    (6) Some codes may need to dump the opacity arrays with the
  59.        frequency as the first index. The freqfirst flag handles this.
  60.  */
  61.  
  62. func multio(filename, opac=, emiss=, srcf=, oscale=, escale=,
  63.         gb=, gav=, gexist=, gscale=, zonelist=, zoneuse=,
  64.         tscale=, noextrap=, freqfirst=)
  65. /* DOCUMENT mf= multio(filename)
  66.          or mf= multio(file)
  67.      opens file FILENAME for use with the multi_streak function.
  68.      The file MUST be subsequently closed using multic, since
  69.      this function produces a hidden reference to the file.  The function
  70.      multif can be used to return an ordinary file pointer, given the
  71.      returned MF structure.  If the argument is already a stream FILE,
  72.      that file will be used.  The call still produces a hidden copy of
  73.      FILE, so you may set your copy of the FILE variable to [], but do
  74.      not close the file.
  75.  
  76.      The following keywords can be used to allow for variations in the
  77.      variable names or units, and to specify the correspondence between
  78.      the zones in this file, and the zones in the master file:
  79.  
  80.      zonelist=index_list
  81.      -or- zonelist=zonelist_name
  82.        is an index list into the (rt,zt) mesh arrays of the master file.
  83.        If ireg is the region number array (having the same dimensions as
  84.        rt or zt, and with its first row and column all 0), and if FILENAME
  85.        contains opacity data only for zones with region numbers 1 and 2,
  86.        you could open the file using:
  87.           mf= multio(filename, zonelist=where(ireg==1 | ireg==2))
  88.        The zonelist should be nil only if the spatial dimensions of the
  89.        opacity and emissivity in this file exactly match those of rt or
  90.        zt in the master file.
  91.        If zonelist is a string, it replaces the default name for the
  92.        zonelist variable stored in the file (see multi_zonelist).
  93.  
  94.      zoneuse=index_list
  95.        The zonelist specifies how the zones in this file correspond
  96.        with those in the master file.  The zoneuse list allows you
  97.        to specify that only some of the zones actually present in the
  98.        opacity and emissivity arrays of this file are to contribute
  99.        to the total.  This might be necessary to avoid double counting
  100.        in a region covered by more than one file.  Hence zoneuse is
  101.        a list of indices into the spatial dimension(s) of the opacity
  102.        and emissivity arrays in this file.  If nil, all zones in this
  103.        file will contribute.  If present, and if zonelist is supplied
  104.        as an array (rather than out of the file), zonelist should
  105.        have the same length as zoneuse.
  106.        As a special case, if zoneuse is a scalar 0, no opacity or
  107.        emissivity will come from this file; this makes sense only if
  108.        this is the master file.
  109.  
  110.      opac=oname, emiss=ename, srcf=sname
  111.        specify non-default names for the opacity, emissivity, and
  112.        osource function arrays.  The defaults are given by the global
  113.        variables mutli_opac, multi_emiss, and multi_srcf (see help).
  114.        If the emissivity array is present in the file, it is preferred
  115.        to the source function array, which will then be ignored.
  116.  
  117.      oscale=opacity_unit, escale=emissivity_unit
  118.        are optional conversion factors to bring the units of the
  119.        opac and emiss (or srcf) arrays into agreement among the various
  120.        files which are to be used in a single run.  The default value
  121.        is 1.0 (i.e.- all files are expected to have the same units).
  122.  
  123.      gb=gbname, gav=gavname, gexist=gexistname
  124.      -or- gexist=group_existence_map
  125.        specify non-default names for the group boundary, group energy,
  126.        and group existence arrays.  The defaults are given by the global
  127.        variables mutli_gb, multi_gav, and multi_gexist (see help).
  128.        If the group boundary array is present in the file, it is preferred
  129.        to the group energy array, which will then be ignored.  The file
  130.        should specify group boundaries if its opacity and emissivity are
  131.        averaged over finite width bins; group energies if its opacity
  132.        and emissivity are computed at points.  The group existence map,
  133.        if present, allows several disjoint spectral regions to exist in
  134.        a single file.  If the data type of gexist is not "string", it
  135.        should be an array of length one less than gb, if gb is present,
  136.        or gav, otherwise.  By this means you can ignore spectral regions
  137.        which are present in the file.
  138.  
  139.      gscale=photon_energy_unit
  140.        is an optional conversion factor to bring the units of the
  141.        gb (or gav) arrays into agreement among the various files
  142.        which are to be used in a single run.  The default value
  143.        is 1.0 (i.e.- all files are expected to have the same units).
  144.  
  145.      tscale=time_unit
  146.        is an optional conversion factor to bring the units of the
  147.        time into agreement among the various files which are to be used
  148.        in a single run.  The default value is 1.0 (i.e.- all files are
  149.        expected to have the same units).
  150.  
  151.      noextrap=1
  152.        if present and non-zero prevents the opacity and emissivity
  153.        data from this file from being extrapolated as 1/hnu^3 in
  154.        master bins at energies above the highest energy bin in this
  155.        file.
  156.  
  157.      freqfirst=0
  158.        if present and non-zero means the frequency index is first
  159.        for the opacity and emissivity arrays, instead of the
  160.        default of frequency index last.
  161.  
  162.    SEE ALSO: multic, multif, multi_streak, MultiFile
  163.              multi_opac, multi_emiss, multi_srcf, multi_gb, multi_gav,
  164.          multi_zonelist
  165.  */
  166. {
  167.   if (!is_stream(filename)) {
  168.     f= openb(filename);
  169.     if (!is_stream(f)) error, "couldn't openb "+filename;
  170.   } else {
  171.     f= filename;
  172.   }
  173.  
  174.   /* create a symbol and stick a hidden reference to the file in it --
  175.      this is as close as Yorick gets to a pointer to a file variable
  176.      :::should re-implement this with lists for Yorick 1.2::: */
  177.   extern _multi_symbols;
  178.   if (is_void(_multi_symbols)) _multi_symbols= array(string, 8);
  179.   list= where(!_multi_symbols);
  180.   if (!numberof(list)) {
  181.     grow, _multi_symbols, array(string, 8);
  182.     list= where(!_multi_symbols);
  183.   }
  184.   i= list(1);
  185.   _multi_symbols(i)= name= "_m_f_"+pr1(i);
  186.   symbol_set, name, f;
  187.  
  188.   mf= MultiFile(f=name);
  189.   if (!is_void(zonelist) && structof(zonelist)!=string) {
  190.     mf.zones= &zonelist;
  191.     zonelist= [];
  192.   }
  193.   if (!is_void(zoneuse)) mf.zuse= &zoneuse;
  194.   if (!is_void(opac)) mf.oname= opac;
  195.   if (!is_void(emiss)) mf.ename= emiss;
  196.   if (!is_void(srcf)) mf.sname= srcf;
  197.   if (!is_void(escale)) mf.escale= escale;
  198.   if (!is_void(oscale)) mf.oscale= oscale;
  199.   if (!is_void(noextrap)) mf.noextrap= noextrap;
  200.   if (!is_void(freqfirst)) {
  201.     mf.freqfirst= freqfirst;
  202.   } else {
  203.     /* default is zero so q-files will work */
  204.     mf.freqfirst= 0;
  205.   }
  206.   _multi_names, mf, gbname=gb, gavname=gav, gexistname=gexist,
  207.                     gscale=gscale, zonesname=zonelist;
  208.  
  209.   if (!is_void(tscale) && tscale>0.0) mf.tscale= tscale;
  210.   else mf.tscale= 1.0;
  211.  
  212.   return mf;
  213. }
  214.  
  215. _multi_symbols= [];
  216.  
  217. func multic(mf)
  218. /* DOCUMENT multic, mf
  219.          or multic, [mf1, mf2, mf3, ...]
  220.      closes a MultiFile created with multif.
  221.      Presented with an array of multifiles, closes them all.
  222.    SEE ALSO: multio, multif
  223.  */
  224. {
  225.   n= numberof(mf);
  226.   for (i=1 ; i<=n ; i++) {
  227.     name= mf(i).f;
  228.     mf(i).f= "?";    /* prevent future references */
  229.     if (!is_void(_multi_symbols)) {
  230.       list= where(_multi_symbols==name);
  231.       if (numberof(list)) _multi_symbols(list)= string(0);
  232.     }
  233.     /* discard hidden file reference */
  234.     symbol_set, name, [];
  235.     /* free associated storage as well */
  236.     mf.zones= mf.gb= mf.gav= mf.gexist= mf.gx= mf.nu= &[];
  237.   }
  238. }
  239.  
  240. func multif(mf)
  241. /* DOCUMENT multif(mf)
  242.      returns an ordinary file pointer for the MultiFile MF.
  243.      Do not use close to close this pointer; just set it to [] when
  244.      you are done.  Use multic to properly close the MF.
  245.    SEE ALSO: multio, multic
  246.  */
  247. {
  248.   return symbol_def(mf.f);
  249. }
  250.  
  251. func multi_streak(mf, rays, slimits, gb=)
  252. /* DOCUMENT result= multi_streak(mf, rays, slimits, gb=common_bins)
  253.      like the streak function, but allows opacity to be built up from
  254.      "slave files", in addition to the "master file" MF(1).  The MF
  255.      parameter is an array of MultiFiles, each created by multif.
  256.      The master file MF(1) contains the mesh, and the master list of dump
  257.      times.  Only dump times which are present in this master list, and
  258.      in every slave file, will be processed.
  259.      The master file MF(1) need not contain any opacity or emissivity data
  260.      at all; each of the slave files MF(2:0) must contain data for at
  261.      least one zone.
  262.  
  263.      The emissivities and opacities from each file are interpolated onto
  264.      a common group structure.  This common group structure can be
  265.      provided via the GB keyword to multi_streak.  If it is not provided,
  266.      GB is computed by examining the group boundary (or center) arrays
  267.      from the master and every slave file, and building a group structure
  268.      which is at least as fine as every component group structure, at every
  269.      point in the spectrum.
  270.  
  271.      Example:
  272.        File family "prob_p00" contains the mesh and opacities and
  273.        emissivities for all zones.  Family "pp_h00" contains post
  274.        processed opacities and emissivities on a much finer spectral
  275.        mesh, but only for zones in regions 1 and 2 of the original
  276.        problem.  You want to transport the emission from the
  277.        inner regions 1 and 2 through the overlying material:
  278.  
  279.          restore, openb("prob_p00"), ireg;
  280.      master= multif("prob_p00", zoneuse=where(ireg>2));
  281.      slave= multif("pp_h00", zonelist=where(ireg==1|ireg==2));
  282.      rays= ...
  283.      slimits= ...
  284.      drat_start= ...
  285.      drat_stop= ...
  286.      result= multi_streak([master,slave], rays, slimits);
  287.      multic, master;
  288.      multic, slave;
  289.  
  290.    SEE ALSO: multio, multic, multif, MultiFile
  291.              multi_opac, multi_emiss, multi_srcf, multi_gb, multi_gav,
  292.          multi_zonelist
  293.          multi_times   (get list of times used by multi_streak)
  294.          multi_bins    (get bins used by multi_streak)
  295.  */
  296. {
  297.   f= multif(mf(1));       /* master file is first one */
  298.  
  299.   /* initialize local variables which are hidden arguments to
  300.      _multi_... functions */
  301.   _multi_files= mf;
  302.   _multi_gb= gb;
  303.   _multi_nchunks= 0;     /* set in _multi__init */
  304.  
  305.   /* set up hook functions for streak */
  306.   streak_times= multi_times;
  307.   _multi_init= _multi__init;
  308.   drat_integrate= _multi_integrate;
  309.   jt= _multi_jt;
  310.  
  311.   return streak(f, rays, slimits);
  312. }
  313.  
  314. func multi_streak_save(outname, mf, rays, slimits, gb=)
  315. /* DOCUMENT multi_streak_save, outname, mf, rays, slimits, gb=common_bins
  316.          or multi_streak_save, outfile, mf, rays, slimits, gb=common_bins
  317.      like the streak function, but allows opacity to be built up from
  318.      "slave files", in addition to the "master file" MF(1) and 
  319.      saves the streak in a PDB history file.  The MF parameter
  320.      is an array of MultiFiles, each created by multif.
  321.      The master file MF(1) contains the mesh, and the master list of dump
  322.      times.  Only dump times which are present in this master list, and
  323.      in every slave file, will be processed.
  324.      The master file MF(1) need not contain any opacity or emissivity data
  325.      at all; each of the slave files MF(2:0) must contain data for at
  326.      least one zone.
  327.  
  328.      If the first argument is OUTFILE, a file variable instead of a
  329.      file name, then that file is used for output.  You can create
  330.      OUTFILE and add static variables to it with save (but do NOT call
  331.      add_record) which streak_save otherwise wouldn't know about.
  332.  
  333.      The output file has history records at the same times as the
  334.      input file.  Each record contains "time" (a double scalar),
  335.      and the two arrays "transp", the transparency (between 0 and 1),
  336.      and "selfem", the self emission (which has the same units as
  337.      ekap in the file F).  The dimensions of transp and selfem
  338.      are ngroup-by-2-by-nrays (where nrays represents zero or more
  339.      dimensions, copied from the RAYS input array).  The RAYS and
  340.      SLIMITS inputs are placed into the output file as non-record
  341.      variables, and any variables in the drat_static option are
  342.      copied form F to the output file.  The gb and gav variables
  343.      are copied from F into the output file as well.  If the drat_glist
  344.      option is present, that is stored in the output file also.
  345.  
  346.      The emissivities and opacities from each file are interpolated onto
  347.      a common group structure.  This common group structure can be
  348.      provided via the GB keyword to multi_streak.  If it is not provided,
  349.      GB is computed by examining the group boundary (or center) arrays
  350.      from the master and every slave file, and building a group structure
  351.      which is at least as fine as every component group structure, at every
  352.      point in the spectrum.
  353.  
  354.      Example:
  355.        File family "prob_p00" contains the mesh and opacities and
  356.        emissivities for all zones.  Family "pp_h00" contains post
  357.        processed opacities and emissivities on a much finer spectral
  358.        mesh, but only for zones in regions 1 and 2 of the original
  359.        problem.  File "prob_strk" contains the streak history. 
  360.        You want to transport the emission from the
  361.        inner regions 1 and 2 through the overlying material:
  362.  
  363.          restore, openb("prob_p00"), ireg;
  364.      master= multif("prob_p00", zoneuse=where(ireg>2));
  365.      slave= multif("pp_h00", zonelist=where(ireg==1|ireg==2));
  366.      fout= openb("prob_strk");
  367.      save, fout, kmax, lmax;
  368.      rays= ...
  369.      slimits= ...
  370.      drat_start= ...
  371.      drat_stop= ...
  372.      result= multi_streak_save(fout, [master,slave], rays, slimits);
  373.      multic, master;
  374.      multic, slave;
  375.  
  376.    SEE ALSO: multio, multic, multif, MultiFile, multi_streak
  377.              multi_opac, multi_emiss, multi_srcf, multi_gb, multi_gav,
  378.          multi_zonelist
  379.          multi_times   (get list of times used by multi_streak)
  380.          multi_bins    (get bins used by multi_streak)
  381.  */
  382. {
  383.   dummy= use_origins(0);
  384.  
  385.   /* create the output file */
  386.   { local f_save; }    /* used in streak_saver
  387.                           -- must NOT be local variable in streak */
  388.   if (is_stream(outname)) f_save= outname;
  389.   else f_save= createb(outname);
  390.  
  391.   f= multif(mf(1));       /* master file is first one */
  392.  
  393.   /* initialize local variables which are hidden arguments to
  394.      _multi_... functions */
  395.   _multi_files= mf;
  396.   _multi_gb= gb;
  397.   _multi_nchunks= 0;     /* set in _multi__init */
  398.  
  399.   /* set up hook functions for streak */
  400.   streak_times= multi_times;
  401.   _multi_init= _multi__init;
  402.   drat_integrate= _multi_integrate;
  403.   jt= _multi_jt;
  404.  
  405.   /* save rays and slimits parameters to output file */
  406.   if (is_void(slimits)) save, f_save, rays;
  407.   else save, f_save, rays, slimits;
  408.  
  409.   /* copy drat_static from input file to output file */
  410.   { extern drat_static; }
  411.   vars= get_vars(f);
  412.   if (drat_gb && is_present(vars, drat_gb)) list= [drat_gb];
  413.   else list= [];
  414.   if (drat_gav && is_present(vars, drat_gav)) grow, list, [drat_gav];
  415.   grow, list, drat_static;
  416.   n= numberof(list);
  417.   for (i=1 ; i<=n ; i++) {
  418.     name= list(i);
  419.     if (!is_present(vars, name))
  420.       error, "\""+name+"\" specified in drat_static option not present";
  421.     value= get_member(f, name);
  422.     add_variable, f_save, -1, name, structof(value), dimsof(value);
  423.     get_member(f_save, name)= value;
  424.   }
  425.  
  426.   /* save drat_glist option if present */
  427.   { extern drat_glist; }
  428.   if (!is_void(drat_glist)) save, f_save, drat_glist;
  429.  
  430.   /* swap out drat_compress and call streak */
  431.     { extern drat_compress; }     /* option */
  432.   drat_compress= streak_saver;
  433.   streak, f, rays, slimits;
  434. }
  435.  
  436. func multi_bins(mf)
  437. /* DOCUMENT multi_bins(mf)
  438.      The MF parameter is an array of MultiFiles, each created by multif.
  439.      Automatically generates the bin structure which will be used by
  440.      multi_streak (if the GB keyword is not specified).
  441.  */
  442. {
  443.   /* get the group structure from the master file, if any */
  444.   mfi= mf(1);
  445.   if (mfi.gb) gb= *mfi.gb;
  446.   else if (mfi.gav) gb= multi_bav(*mfi.gav);
  447.  
  448.   /* build an appropriate bin boundary array */
  449.   n= numberof(mf);
  450.   for (i=2 ; i<=n ; i++) {
  451.     mfi= mf(i);
  452.     if (mfi.gb) gbf= *mfi.gb;
  453.     else if (mfi.gav) gbf= multi_bav(*mfi.gav);
  454.     else error, "no gb or gav arrays in slave file "+pr1(i);
  455.     if (is_void(gb)) gb= gbf;
  456.     else gb= _multi_bins(0, gb, gbf);
  457.   }
  458.  
  459.   return gb;
  460. }
  461.  
  462. func _multi_jt(f, time)
  463. {
  464.   _jt, f, time;
  465.   n= numberof(_multi_files);
  466.   for (i=2 ; i<=n ; i++) {
  467.     mfi= _multi_files(i);
  468.     _jt, multif(mfi), time/mfi.tscale;
  469.   }
  470. }
  471.  
  472. /* streak_times is "swapped out" in multi_streak */
  473. _orig_streak_times= streak_times;
  474.  
  475. func multi_times(f)
  476. /* DOCUMENT times= multi_times(mf)
  477.      returns the list of times which will be used by multi_streak.  This
  478.      is the subset of streak_times(mf(1)) which occur in all of the slave
  479.      files.  The drat_start and drat_stop times work as usual.
  480.  */
  481. {
  482.   if (!is_stream(f)) {
  483.     mf= f;
  484.     f= multif(mf(1));
  485.   } else {
  486.     mf= _multi_files;  /* if called as streak_times from streak function,
  487.               _multi_files is local to multi_streak */
  488.   }
  489.   if (!is_void(mf) &&
  490.       abs(mf(1).tscale-1.0)>1.e-6) error, "tscale illegal in master file";
  491.   times= _orig_streak_times(f);
  492.  
  493.   n= numberof(mf);
  494.   for (i=2 ; i<=n ; i++) {
  495.     mfi= mf(i);
  496.     times= _multi_times(times, mfi.tscale*get_times(multif(mfi)));
  497.   }
  498.  
  499.   /* give the multi package a chance to initialize just before the main
  500.      loop over time is entered in the streak function */
  501.   if (!is_void(_multi_init)) _multi_init;
  502.  
  503.   return times;
  504. }
  505.  
  506. local multi_memory;
  507. /* DOCUMENT multi_memory
  508.      amount of memory used to determine size of spectral chunks.
  509.      Default is 2000000, which keeps the memory required per chunk
  510.      to under a few megabytes.
  511.  */
  512. multi_memory= 2000000;
  513.  
  514. func _multi__init
  515. /* xxDOCUMENT multi_init
  516.      is called immediately before the streak function enters its loop
  517.      over time.  It is responsible for computing all of the photon group
  518.      interpolation weights for the master and slave files.
  519.  */
  520. {
  521.   extern _multi_files, _multi_gb;  /* local to multi_streak */
  522.  
  523.   /* build an appropriate bin boundary array if necessary */
  524.   if (is_void(_multi_gb)) _multi_gb= multi_bins(_multi_files);
  525.  
  526.   /* compute a reasonable number of chunks, trying to limit memory
  527.      usage to a few megabytes */
  528.   extern _multi_nchunks;    /* local to multi_streak */
  529.   kmx= dimsof(get_member(multif(_multi_files(1)), drat_rt));
  530.   lmx= kmx(3);
  531.   kmx= kmx(2);
  532.   ngroup= numberof(_multi_gb)-1;
  533.   _multi_nchunks= min(max(multi_memory/(8*8*kmx*lmx), 1), ngroup);
  534.   if (ngroup%_multi_nchunks) _multi_nchunks= ngroup/_multi_nchunks + 1;
  535.   else _multi_nchunks= ngroup/_multi_nchunks;
  536.  
  537.   /* compute the various indices and other data required for the
  538.      interpolation and splitting into chunks */
  539.   n= numberof(_multi_files);
  540.   for (i=1 ; i<=n ; i++)
  541.     _multi_chunk, i, _multi_gb, _multi_nchunks;
  542. }
  543.  
  544. func _multi_integrate(f, mesh, time, irays, slimits)
  545. /* DOCUMENT atten_emit= _multi_integrate(f, mesh, time, irays, slimits)
  546.      is the default drat_integrate routine.
  547.      On entry, file F is positioned at TIME, from which MESH has already
  548.      been read.  IRAYS and SLIMITS are the rays coordinates (in internal
  549.      format) and integration limits.
  550.      The result should be ngroup-by-2-by-raydims, where the second index
  551.      is 1 for the attenuation factor, 2 for the self-emission (specific
  552.      intensity due to emission along the ray).
  553.    OPTIONS: drat_linear, drat_ocompute, drat_oadjust,
  554.             drat_emult, drat_amult, drat_omult, drat_nomilne,
  555.         drat_ekap, drat_akap, drat_glist
  556.    SEE ALSO: streak, multi_streak
  557.  */
  558. {
  559.   /* build arrays which will hold the opacity and emissivity chunks */
  560.   ngroup= numberof(_multi_gb)-1;
  561.   size= ngroup/_multi_nchunks;
  562.   nlong= ngroup%_multi_nchunks;
  563.   if (nlong) size++;
  564.   /* local here, external to _multi_accum */
  565.   opac= emiss= array(0.0, dimsof(get_member(f,drat_rt)), size);
  566.  
  567.   /* select the appropriate integrator */
  568.   integrator= drat_linear? integ_linear : integ_flat;
  569.  
  570.   /* loop over spectral chunks */
  571.   nfiles= numberof(_multi_files);
  572.   n2= 0;
  573.   for (i=1 ; i<=_multi_nchunks ; i++) {
  574.     n1= n2+1;
  575.     n2+= size;
  576.  
  577.     /* accumulate opacities and emissivities from all files */
  578.     for (j=1 ; j<=nfiles ; j++)
  579.       _multi_accum, _multi_files(j), _multi_gb(n1:n2+1), i;
  580.  
  581.     /* do transport integral for these rays cutting through current mesh */
  582.     emiss/= opac+1.e-99;  /* do it here to reduce number of temporaries */
  583.  
  584.     /* apply absorption and emission multipliers */
  585.     { extern drat_amult, drat_emult, drat_omult; }   /* options */
  586.     if (!is_void(drat_emult)) emiss*= drat_emult;
  587.     if (!is_void(drat_amult)) {
  588.       opac*= drat_amult+1.e-20;
  589.       emiss/= drat_amult+1.e-20;
  590.     }
  591.     if (!is_void(drat_omult)) opac*= drat_omult;
  592.  
  593.     /* point center the source function (in place) if requested */
  594.     { extern drat_linear, drat_nomilne; }   /* options */
  595.     if (drat_linear) {
  596.       if (structof(emiss)!=double) emiss= double(emiss);
  597.       pcen_source, opac, emiss, mesh, drat_nomilne;
  598.     }
  599.  
  600.     chunk= integrator(opac, emiss, irays, mesh, slimits);
  601.  
  602.     if (i==1) {
  603.       /* result is same as chunk except in spectral dimension,
  604.      which is length ngroup instead of chunk size */
  605.       dims= dimsof(chunk);
  606.       dims(2)= ngroup;
  607.       instant= array(0.0, dims);
  608.     }
  609.     instant(n1:n2,..)= chunk;
  610.     chunk= [];  /* free memory ASAP */
  611.  
  612.     if (i==nlong) {
  613.       /* ngroup%_multi_nchunks chunks are length ngroup/_multi_nchunks+1,
  614.      the rest are length ngroup/_multi_nchunks
  615.      Adjust size for the next pass -- note that nlong<_multi_nchunks
  616.      by construction, so this can never happen on the last pass; it
  617.      never happens at all if nlong==0.  */
  618.       size--;
  619.       dims= dimsof(opac);
  620.       dims(0)= size;
  621.       opac= emiss= [];  /* destroy before recreating */
  622.       opac= emiss= array(0.0, dims);
  623.     } else {
  624.       opac(..)= emiss(..)= 0.0;  /* initialize for next pass */
  625.     }
  626.   }
  627.  
  628.   return instant;
  629. }
  630.  
  631. /* ------------------------------------------------------------------------ */
  632.  
  633. struct MultiFile {
  634.   string f;         /* name of the file variable for this file
  635.                -- no such thing as address of a file variable */
  636.  
  637.   string oname;     /* name of the opacity variable in f
  638.                -- default is "opac" if present in the file */
  639.   string ename;     /* name of the emissivity variable in f, or 0 if
  640.                source function used instead
  641.                -- default is "emiss" if present in the file */
  642.   string sname;     /* name of the source function variable in f,
  643.                ignored if ename!=0
  644.                -- default is "srcf" if present in the file */
  645.  
  646.   double oscale;    /* scale factor for correcting opacity units */
  647.   double escale;    /* scale factor for correcting emissivity units
  648.                -- will be used to correct srcf*opac if srcf used */
  649.  
  650.   pointer zuse;     /* pointer to list of zone indices to use in this file
  651.                -- If &[], all zones in this file are used.
  652.                -- As a special case, if &array(0), use no zones.  */
  653.  
  654.   pointer zones;    /* list of master file zone indices for the zones which
  655.                are present in this post-processing file
  656.                -- The length of this list must be the same as the
  657.                product of the non-spectral dimensions of the opacity
  658.                and emissivity arrays, or the same length as *zuse,
  659.                if that array is present.
  660.                -- zones==0 means that all zones in the master file
  661.                (including zeroes in the first row and column) are
  662.                present in this file.  */
  663.  
  664.   /* This standard interface supports two spectral models:
  665.      (1) Zone-centered model: Opacities and emissivities are defined as
  666.          integrals over photon energy bins, so that the group structure
  667.      is specified by a bin boundary energy list gb.
  668.      (2) Point-centered model: Opacities and emissivities are defined at
  669.          particular discrete photon energies, so that the group structure
  670.      is specified by a photon energy list gav.
  671.      In case (1), gav==0, and in case (2), gb==0.
  672.      In either case, gb or gav must be a strictly increasing list of
  673.      photon energies.  However, a single file may contain data in several
  674.      disjoint spectral regions.  This information, together with spectral
  675.      chunking information, is stored in the index array, nu.
  676.    */
  677.   pointer gb;       /* bin boundary array list for this file */
  678.   pointer gav;      /* photon energy array list for this file */
  679.   pointer gexist;   /* group existence map (if gb!=0) or connection map
  680.                (if gav!=0), or 0 if all groups exist or all
  681.                photon energies are connected */
  682.  
  683.   pointer gx;
  684.   /* (*gx) is a pointer to a n_spectral_regions -by- 2 array
  685.      of extreme values of each spectral region, used to limit the
  686.      interpolation boundaries.  (*gx)(,1) is the lower boundary,
  687.      while (*gx)(,2) is the upper boundary.  */
  688.  
  689.   pointer nu;
  690.   /* (*nu) is 6 -by- n_spectral_regions -by- n_chunks index array:
  691.      (*nu)(, spectral-region, chunk)
  692.      can be written [i1,i2,i3,i4,i5,i6] where:
  693.  
  694.      gb(i1:i2+1) or gav(i1:i2) are the bins required to compute this
  695.           chunk in this spectral region.
  696.       i1==0 if, and only if, this file make no contribution at all
  697.       to this chunk in this spectral region.
  698.      multi_gb(i3:i4) are the master bins affected by this chunk
  699.      If i6!=0, opacity(i4+1:i5) will be extrapolated as 1/nu^3 from the
  700.      value at gb(i6:i6+1) or gav(i6).  This will be done even if i1==0.  */
  701.  
  702.   double tscale;    /* scale factor for correcting time units */
  703.  
  704.   int noextrap;  /* 1/nu^3 extrapolation will be done unless non-zero */
  705.  
  706.   int freqfirst;  /* frequency index is last unless non-zero */
  707. }
  708.  
  709. func _multi_chunk(i, gb, nchunks)
  710. /* xxDOCUMENT _multi_chunk, i, gb, nchunks
  711.      computes the nu and gi members of the MultiFile struct _multi_files(I),
  712.      so that the master bin boundary array GB can be computed in NCHUNKS
  713.      chunks.
  714.  */
  715. {
  716.   mf= _multi_files(i);
  717.   /* get group structure out of file */
  718.   gbf= *mf.gb;
  719.   gavf= *mf.gav;
  720.   gexist= *mf.gexist;
  721.   if (is_void(gbf) && is_void(gavf)) return;
  722.  
  723.   /* compute chunk boundaries */
  724.   nb= numberof(gb);
  725.   ngroup= nb-1;
  726.   if (nchunks>ngroup) nchunks= ngroup;
  727.   chunk= array(ngroup/nchunks, nchunks);
  728.   n= ngroup%nchunks;
  729.   if (n) chunk(1:n)+= 1;
  730.   chunk= chunk(cum);
  731.   cbot= chunk(1:-1)+1;
  732.   ctop= chunk(2:0);
  733.  
  734.   /* get j1= indices of beginnings of spectral regions,
  735.          j2= indices of ends of spectral regions
  736.      in file bins */
  737.   j1= _multi_spectreg(gexist, gbf, gavf);
  738.   j2= j1(,2);
  739.   j1= j1(,1);
  740.  
  741.   /* get master bin numbers containing j1 and j2
  742.      -- after this, bot and top are bin number containing the
  743.         bottom and top of each of the spectral regions
  744.      -- 1 and ngroup+2 are still possible if the region partially
  745.         hangs over the master bin structure */
  746.   if (!is_void(gbf)) {
  747.     gf= gbf;
  748.     nf= max(j2)-1;   /* top group index in file */
  749.   } else {
  750.     gf= gavf;
  751.     nf= max(j2);
  752.   }
  753.   gbot= gf(j1);
  754.   gtop= gf(j2+(!is_void(gbf)));
  755.   bot= digitize(gbot, gb)-1;
  756.   top= digitize(-gtop, -gb)-1;  /* if one of the gtop is exactly equal to
  757.                    one of the gb, want the smaller bin
  758.                    for top (see digitize help) */
  759.  
  760.   /* get the master group indices of the minimum i3 and maximum i4 group
  761.      which can be affected by each spectral region for each chunk
  762.      -- note that i3>=1 and i4<=ngroup are guaranteed after this, so
  763.         gb(i3) and gb(i4+1) below will never blow up */
  764.   i3= max(bot, cbot(-,));
  765.   i4= min(top, ctop(-,));
  766.  
  767.   /* get minimum i1 and maximum i2 bin indices in file
  768.      which must be extracted for each spectral region and each chunk
  769.      -- note that i3>i4 in empty spectral regions, which will be
  770.         removed later */
  771.   i1= digitize(gb(i3), gf)-1;
  772.   i2= digitize(-gb(i4+1), -gf);  /* see digitize comment above */
  773.   if (!is_void(gbf)) i2--;  /* top group one less than top bin boundary */
  774.   /* can't retrieve anything beyond min(gf) and max(gf) */
  775.   i1= max(1, i1);
  776.   i2= min(nf, i2);
  777.  
  778.   /* zero empty spectral region/chunks -- they haven't been computed
  779.      correctly anyway */
  780.   list= where(i3>i4);
  781.   if (numberof(list)) i1(list)= i2(list)= i3(list)= i4(list)= 0;
  782.  
  783.   /* find spectral regions where extrapolation is necessary */
  784.   cbot-= 1;
  785.   cbot= cbot(-:1:numberof(j2),);
  786.   i5= i6= 0*i1;
  787.   list= where(top<ctop(-,));  /* see definition of i4 above */
  788.   if (!mf.noextrap && numberof(list)) {
  789.     i6(list)= j2(,-:1:numberof(ctop))(list);
  790.     i5(list)= ctop(-:1:numberof(j2),)(list);
  791.     i4e= i4(list);
  792.     l= where(i4e==0);
  793.     if (numberof(l)) {
  794.       /* fix any chunks which require extrapolation, but which are
  795.      completely above the spectral region being extrapolated */
  796.       i4e(l)= cbot(list)(l);
  797.       i4(list)= i4e;
  798.     }
  799.   }
  800.  
  801.   /* adjust i3, i4, and i5 to apply to the current chunk, rather
  802.      than to the entire gb array */
  803.   list= where(i3);
  804.   if (numberof(list)) i3(list)-= cbot(list);
  805.   list= where(i4);
  806.   if (numberof(list)) i4(list)-= cbot(list);
  807.   list= where(i5);
  808.   if (numberof(list)) i5(list)-= cbot(list);
  809.  
  810.   /* install index lists into mf */
  811.   _multi_files(i).nu= &transpose([i1,i2,i3,i4,i5,i6],2);
  812.  
  813.   /* generate the interpolation bin limits for each spectral region
  814.      and install into mf */
  815.   _multi_files(i).gx= &[gbot, gtop];
  816. }
  817.  
  818. func _multi_spectreg(gexist, gbf, gavf)
  819. /* xxDOCUMENT _multi_spectreg(gexist, gbf, gavf)
  820.      crack GBF (if non-nil) or GAVF (if GBF nil) into spectral regions
  821.      according to GEXIST.  Returns 2-by-n_spectral_regions array of
  822.      minimum and maximum group indices for each spectral region.
  823.      In the case of GBF, GEXIST is a bin existence map, dimension one
  824.      less than GBF, non-zero in bins which exist.
  825.      In the case of GAVF, GEXIST is a photon connectivity array, dimension
  826.      one less than GAVF, zero between consecutive elements of GAVF which
  827.      are not connected.  Single, isolated points are removed entirely.
  828.      In either case, GAVF or GBF must be a strictly increasing list.
  829.  */
  830. {
  831.   n= numberof(gbf);
  832.   if (is_void(gexist)) {
  833.     j1= [1];
  834.     j2= [(n? n-1 : numberof(gavf))];
  835.   } else {
  836.     gexist= grow([0], (gexist!=0), [0]);
  837.     gexist= gexist(dif);
  838.     j1= where(gexist>0);  /* off to on transitions */
  839.     j2= where(gexist<0);  /* on to off transitions */
  840.     if (numberof(j2)) {
  841.       if (n) {
  842.     j2--;   /* bin index, not boundary index */
  843.       } else if (anyof(j2==j1+1)) {
  844.     list= where(j2>j1+1);
  845.     j1= j1(list);
  846.     j2= j2(list);
  847.       }
  848.     }
  849.     if (numberof(j2)<1) error, "group existence map for file is void";
  850.   }
  851.   return [j1,j2];
  852. }
  853.  
  854. func _multi_accum(mf, gb, chunk)
  855. /* xxDOCUMENT _multi_accum, mf, gb, chunk
  856.      increments externals opac, emiss from MultiFile MF, chunk number
  857.      CHUNK.  GB is the master bin boundary array.
  858.  */
  859. {
  860.   /* unpack the mf MultiFile struct */
  861.   nu= mf.nu;
  862.   if (!nu) return;         /* detect no-op quickly */
  863.   nu= (*nu)(,,chunk);
  864.   if (noneof(nu)) return;  /* detect no-op quickly */
  865.   f= multif(mf);
  866.   oname= mf.oname;
  867.   ename= mf.ename;
  868.   sname= mf.sname;
  869.   oscale= mf.oscale;
  870.   escale= mf.escale;
  871.   freqfirst= mf.freqfirst;
  872.   zuse= *mf.zuse;   /* if scalar 0, nu==&[] above, and won't get here */
  873.   zones= *mf.zones;
  874.   gbf= *mf.gb;
  875.   gavf= *mf.gav;
  876.   gx= (*mf.gx);
  877.   gbot= gx(,1);
  878.   gtop= gx(,2);
  879.  
  880.   n= dimsof(nu)(3);
  881.   for (i=1 ; i<=n ; i++) {
  882.     /* loop over disjoint spectral regions in this file
  883.        -- this is a little inefficient, as presumably at most two
  884.           adjacent ones will have non-zero i1 */
  885.     i1= nu(1,i);
  886.     i2= nu(2,i);
  887.     i3= nu(3,i);
  888.     i4= nu(4,i);
  889.     i5= nu(5,i);
  890.     i6= nu(6,i);
  891.  
  892.     if (i1) {
  893.       /* compute integral of opacity and emissivity over this chunk,
  894.          interpolated onto the master bin structure */
  895.       gbsub= gb(i3:i4+1);
  896.       gbm= min(max(gbsub,gbot(i)),gtop(i));
  897.       /* NOTE: oraw is two-dimensional */
  898.       if(freqfirst) {
  899.     oraw= transpose(get_member(f,oname)(i1:i2,*));
  900.       } else {
  901.     oraw= get_member(f,oname)(*,i1:i2);
  902.       }
  903.       if (!is_void(gbf)) {
  904.     /* zone-centered spectral model */
  905.     g12= gbf(i1:i2+1);
  906.     dg12= g12(-,dif);
  907.     if (is_void(gavf)) o= interp((oraw*dg12)(,cum), g12, gbm, 2);
  908.     else o= integ(oraw, gavf(i1:i2), gbm, 2);
  909.     if (ename) {
  910.       oraw= [];
  911.       if(freqfirst) {
  912.         e= transpose(get_member(f,ename)(i1:i2,*));
  913.         e= interp((e*dg12)(,cum), g12, gbm, 2);
  914.       } else {
  915.         e= interp((get_member(f,ename)(*,i1:i2)*dg12)(,cum), g12, gbm, 2);
  916.       }
  917.     } else {
  918.       if(freqfirst) {
  919.         e= transpose(oraw*get_member(f,sname)(i1:i2,*));
  920.         e= interp((e*dg12)(,cum), g12, gbm, 2);
  921.       } else {
  922.         e= interp((oraw*get_member(f,sname)(*,i1:i2)*dg12)(,cum),
  923.               g12, gbm, 2);
  924.       }
  925.       oraw= [];
  926.     }
  927.       } else {
  928.     /* point-centered spectral model */
  929.     g12= gavf(i1:i2);
  930.     dg12= g12(-,dif);
  931.     o= integ(oraw, g12, gbm, 2);
  932.     if (ename) {
  933.       oraw= [];
  934.       if(freqfirst) {
  935.         e= transpose(get_member(f,ename)(i1:i2,*));
  936.         e= integ(e, g12, gbm, 2);
  937.       } else {
  938.         e= integ(get_member(f,ename)(*,i1:i2), g12, gbm, 2);
  939.       }
  940.     } else {
  941.       if(freqfirst) {
  942.         e= transpose(oraw*get_member(f,sname)(i1:i2,*));
  943.         e= integ(e, g12, gbm, 2);
  944.       } else {
  945.         e= integ(oraw*get_member(f,sname)(*,i1:i2), g12, gbm, 2);
  946.       }
  947.       oraw= [];
  948.     }
  949.       }
  950.  
  951.       /* differentiate and scale */
  952.       dgb= gbsub(-,dif);
  953.       o= o(,dif)/(dgb/oscale);
  954.       e= e(,dif)/(dgb/escale);
  955.  
  956.       /* accumulate into opac and emiss externals */
  957.       if (!is_void(zuse)) {
  958.     /* only a subset of the zones present in this file contribute */
  959.     o= o(zuse,);
  960.     e= e(zuse,);
  961.       }
  962.       if (is_void(zones)) {
  963.     /* all zones represented in this file */
  964.     opac(*,i3:i4)+= o;
  965.     emiss(*,i3:i4)+= e;
  966.       } else {
  967.     /* subset of zones in this file */
  968.     opac(zones,1,i3:i4)+= o;
  969.     emiss(zones,1,i3:i4)+= e;
  970.       }
  971.     }
  972.  
  973.     if (i6) {
  974.       /* must extrapolate opacity as 1/hnu^3 and emissivity as
  975.          exp(-hnu/kT) */
  976.       /* first, get gb3, proportional to the integral of 1/hnu^3 across
  977.      the master bins gb (assumed external) */
  978.       gb3= gb(i4+1:i5+1);
  979.       gb3= gb3(zcen)/(gb3(1:-1)*gb3(2:0))^2;
  980.       if (!is_void(gbf)) gf= avg(gbf(i6:i6+1));
  981.       else gf= gavf(i6);
  982.       gb3*= gf^3;
  983.  
  984.       /* next, get the opacity and emissivity in bin i6 of this file */
  985.       if(freqfirst) {
  986.     o= transpose(get_member(f,oname)(i6,*));
  987.       } else {
  988.     o= get_member(f,oname)(*,i6);
  989.       }
  990.       if (ename) {
  991.     if(freqfirst) {
  992.       e= transpose(get_member(f,ename)(i6,*));
  993.     } else {
  994.       e= get_member(f,ename)(*,i6);
  995.     }
  996.       } else {
  997.     if(freqfirst) {
  998.       e= transpose(o*get_member(f,sname)(i6,*));
  999.     } else {
  1000.       e= o*get_member(f,sname)(*,i6);
  1001.     }
  1002.       }
  1003.  
  1004.       if (!is_void(zuse)) {
  1005.     /* only a subset of the zones present in this file contribute */
  1006.     o= o(zuse);
  1007.     e= e(zuse);
  1008.       }
  1009.       if (is_void(zones)) {
  1010.     /* all zones represented in this file */
  1011.     opac(*,i4+1:i5)+= o*(gb3(-,)*oscale);
  1012.     emiss(*,i4+1:i5)+= e*(gb3(-,)*escale);
  1013.       } else {
  1014.     /* subset of zones in this file */
  1015.     opac(zones,1,i4+1:i5)+= o*(gb3(-,)*oscale);
  1016.     emiss(zones,1,i4+1:i5)+= e*(gb3(-,)*escale);
  1017.       }
  1018.     }
  1019.  
  1020.     /* free temporary storage before next pass */
  1021.     o= e= [];
  1022.  
  1023.   }  /* end of loop over spectral regions */
  1024. }
  1025.  
  1026. local multi_opac, multi_emiss, multi_srcf;
  1027. /* DOCUMENT multi_opac, multi_emiss, multi_srcf
  1028.      are the default names of the opacity, emissivity, and source function
  1029.      arrays in post-processing files.  By default, they are "opac",
  1030.      "emiss", and "srcf", respectively.  If none of these are present,
  1031.      drat_akap and drat_ekap will also be tried for the opacity and
  1032.      source function, respectively.  If emissivity is present, source
  1033.      function will be ignored.
  1034.      The units of opacity are inverse length, of emissivity power per
  1035.      photon energy per sterradian per volume, and source function
  1036.      power per photon energy per sterradian per area.
  1037.  */
  1038.  
  1039. local multi_gb, multi_gav, multi_gexist;
  1040. /* DOCUMENT multi_gb, multi_gav, multi_gexist
  1041.      are the default names of the group boundary, group center, and group
  1042.      existence arrays in post-processing files.  By default, they are
  1043.      "gb", "gav", and "gexist".  If neither is present, drat_gb and
  1044.      drat_gav will also be tried.
  1045.      If present, gb is the photon bin boundary array, and gexist (if
  1046.      present) has one fewer element and is non-zero in bins which exist.
  1047.      Otherwise, gav is the photon energy array, and gexist (if present)
  1048.      has one fewer element and is zero between elements of gav which are
  1049.      not connected.  Isolated points in gav are removed entirely.
  1050.      Either gb or gav must be strictly increasing, and has units of
  1051.      photon energy.
  1052.  */
  1053.  
  1054. local multi_zonelist;
  1055. /* DOCUMENT multi_zonelist
  1056.      is the default name of the variable which is a list of 1-origin
  1057.      zone indices in the mesh of the master file.
  1058.  */
  1059.  
  1060. func _multi_names(mf, gbname=, gavname=, gexistname=, gscale=, zonesname=)
  1061. /* xxDOCUMENT _multi_names, mf, gbname=..., gavname=...,
  1062.                             gexistname=..., gscale=..., zonesname=...
  1063.      find all relevant variable names in the file specified by MF.f.
  1064.      Reads in the gb (or gav) and gexist arrays, setting MF.gb, MF.gav,
  1065.      and MF.gexist appropriately, multiplying the values read from the
  1066.      file by GSCALE if supplied.
  1067.      Reads in MF.zones if not already set and if zonelist is present in
  1068.      the file.
  1069.      If MF.oscale or MF.escale is 0.0, it is set to 1.0; if less than
  1070.      zero, it is set to 0.0.
  1071.      If the routine fails, mf.oname will be string(0) on output.
  1072.  */
  1073. {
  1074.   f= multif(mf);
  1075.   vars= get_vars(f);
  1076.  
  1077.   zuse= *mf.zuse;
  1078.   contrib= (is_void(zuse) || zuse(1));
  1079.  
  1080.   if (contrib) {
  1081.  
  1082.     if (!mf.zones) {
  1083.       if (is_void(zonesname)) zonesname= multi_zonelist;
  1084.       if (is_present(vars, zonesname)) zones= get_member(f, zonesname);
  1085.       if (!is_void(zones) && !is_void(zuse)) zones= zones(zuse);
  1086.       mf.zones= &zones;
  1087.     }
  1088.  
  1089.     if (mf.oname && !is_present(vars, mf.oname)) mf.oname= string(0);
  1090.     if (!mf.oname) {
  1091.       if (is_present(vars, multi_opac)) mf.oname= multi_opac;
  1092.       else if (is_present(vars, drat_akap)) mf.oname= drat_akap;
  1093.     }
  1094.  
  1095.     if (mf.ename && !is_present(vars, mf.ename)) mf.ename= string(0);
  1096.     if (mf.sname && !is_present(vars, mf.sname)) mf.sname= string(0);
  1097.     if (!mf.ename && !mf.sname) {
  1098.       if (is_present(vars, multi_emiss)) mf.ename= multi_emiss;
  1099.       else if (is_present(vars, multi_srcf)) mf.sname= multi_srcf;
  1100.       else if (is_present(vars, drat_ekap)) mf.sname= drat_ekap;
  1101.       else mf.oname= string(0);
  1102.     }
  1103.  
  1104.     if (!mf.gb && !mf.gav && mf.oname) {
  1105.       if (is_void(gbname) && is_void(gavname)) {
  1106.     if (is_present(vars, multi_gb)) gbname= multi_gb;
  1107.     else if (is_present(vars, drat_gb)) gbname= drat_gb;
  1108.     if (is_present(vars, multi_gav)) gavname= multi_gav;
  1109.     else if (is_present(vars, drat_gav)) gavname= drat_gav;
  1110.       } else {
  1111.     if (!is_void(gbname) && !is_present(vars, gbname))
  1112.       error, "gbname= "+gbname+" is not present in file";
  1113.     if (!is_void(gavname) && !is_present(vars, gavname))
  1114.       error, "gavname= "+gavname+" is not present in file";
  1115.       }
  1116.       if (is_void(gexistname)) {
  1117.     if (is_present(vars, multi_gexist)) gexistname= multi_gexist;
  1118.       } else {
  1119.     if (!is_void(gexistname) && !is_present(vars, gexistname))
  1120.       error, "gexistname= "+gexistname+" is not present in file";
  1121.       }
  1122.       if (!gscale) gscale= 1.0;
  1123.  
  1124.       if (!is_void(gexistname)) mf.gexist= &get_member(f, gexistname);
  1125.       if (!is_void(gbname)) mf.gb= &(gscale*get_member(f, gbname));
  1126.       if (!is_void(gavname)) mf.gav= &(gscale*get_member(f, gavname));
  1127.       if (!mf.gb && !mf.gav) mf.oname= mf.ename= string(0);
  1128.     }
  1129.  
  1130.     if (mf.oscale==0.0) mf.oscale= 1.0;
  1131.     else if (mf.oscale<0.0) mf.oscale= 0.0;
  1132.     if (mf.escale==0.0) mf.escale= 1.0;
  1133.     else if (mf.escale<0.0) mf.escale= 0.0;
  1134.  
  1135.   } else {
  1136.     mf.oname= string(0);
  1137.     mf.gb= mf.gav= &[];
  1138.   }
  1139. }
  1140.  
  1141. /* ------------------------------------------------------------------------ */
  1142. /* procedures to assist in combining different sets of photon bins, or
  1143.    otherwise generating photon group structure */
  1144.  
  1145. func multi_bav(gav)
  1146. /* DOCUMENT multi_bav(gav)
  1147.      returns bin boundaries for the bin centers gav.
  1148.      The bin boundaries are taken at the geometric means between
  1149.      consecutive gav(i), with the endpoints extended slightly beyond
  1150.      the endpoints of gav.
  1151.  */
  1152. {
  1153.   gb= log(gav + 1.e-99)(pcen);
  1154.   gb(1)-= gb(2)-gb(1);
  1155.   gb(0)+= gb(0)-gb(-1);
  1156.   return exp(gb);
  1157. }
  1158.  
  1159. func _multi_bins(nfinal, ..)
  1160. /* DOCUMENT gb= _multi_bins(nfinal, gb1, gb2, ...)
  1161.      returns NFINAL+1 boundaries of NFINAL bins constructed by combining
  1162.      the input bin structures GB1, GB2, etc.
  1163.      Use NFINAL=0 to get at least the resolution in the finest GBi in
  1164.      every region of the spectrum.
  1165.  
  1166.      This is done by constructing a total bin density function
  1167.      (#bins/energy width), as the maximum of the bin density of each
  1168.      component.  This total bin density function is integrated, and
  1169.      the integral is divided into NFINAL equal parts; the points in
  1170.      energy at which this division must be made are the returned bin
  1171.      boundaries.
  1172.  
  1173.      In the returned bin structure, the density of bins is everywhere
  1174.      proportional to the densest bins in any of the GBi.
  1175.  */
  1176. {
  1177.   boundary= density= [];
  1178.   while (more_args()) {
  1179.     /* get next input boundary array */
  1180.     gb= next_arg();
  1181.     if (is_void(gb)) continue;
  1182.     gd= 1./gb(dif);  /* bin density is zone centered -- constant
  1183.             between the bin boundaries */
  1184.     nb= numberof(boundary);
  1185.  
  1186.     if (!nb) {
  1187.       /* first pass */
  1188.       boundary= gb;
  1189.       density= gd;
  1190.  
  1191.     } else {
  1192.       /* subsequent passes */
  1193.       /* combine lists of bin boundaries, sort into increasing order,
  1194.      and eliminate any duplicates */
  1195.       combine= grow(boundary, gb);
  1196.       combine= combine(sort(combine));
  1197.       combine= multi_no_dups(combine);
  1198.  
  1199.       /* get lists of indices of the center of each of the new bins in
  1200.      the old boundary and gb arrays */
  1201.       comcen= combine(zcen);
  1202.       list1= digitize(comcen, boundary);
  1203.       list2= digitize(comcen, gb);
  1204.  
  1205.       /* add guard zones for points beyond boundary or gb separately */
  1206.       density= grow([0.], density, [0.]);
  1207.       gd= grow([0.], gd, [0.]);
  1208.  
  1209.       /* the new density is the maximum of the two component densities,
  1210.      and the new boundary is the combination of the two component
  1211.      boundaries */
  1212.       density= max(density(list1), gd(list2));
  1213.       boundary= combine;
  1214.     }
  1215.   }
  1216.   if (is_void(boundary)) return [];
  1217.  
  1218.   /* now have a histogram density(boundary) of the maximum bin density
  1219.      between every energy mentioned in any of the gb -- the integral of
  1220.      this histogram is the piecewise linear function total(boundary): */
  1221.   total= (density*boundary(dif))(cum);
  1222.  
  1223.   /* note that the last point, total(0) is the "natural" number of bins
  1224.      required to get at least the spectral resolution in any of the
  1225.      input gb */
  1226.   if (nfinal<1) nfinal= long(total(0)+0.99);
  1227.  
  1228.   /* should this be enhanced to allow the caller to specify endpoints
  1229.      other than boundary(1) and boundary(0)? */
  1230.   return interp(boundary, total/total(0), span(0.,1.,nfinal+1));
  1231. }
  1232.  
  1233. func multi_line(nbins, hnu0, dhnu, dhnu_min)
  1234. /* DOCUMENT gb= multi_line(nbins, hnu0, dhnu, dhnu_min)
  1235.      returns 2*NBINS+1 bin boundary energies for 2*NBINS bins
  1236.      cenetered around a spectral line at HNU0 of width DHNU.  The
  1237.      result begins at HNU0-DHNU and ends at HNU0+DHNU.  The finest
  1238.      two bins (nearest HNU0) has width DHNU_MIN, and the remaining
  1239.      bins have equal ratio widths as you move away from HNU0.
  1240.  */
  1241. {
  1242.   dhnu= double(dhnu);
  1243.   r= 0.;
  1244.   if (dhnu_min && dhnu_min>0. && dhnu_min<dhnu) {
  1245.     require, "series.i";
  1246.     r= gseries_r(dhnu/dhnu_min, nbins);
  1247.   }
  1248.   if (r>1.) widths= dhnu_min*exp(log(r)*indgen(0:nbins-1));
  1249.   else widths= dhnu/nbins;
  1250.   widths= grow(widths(0:1:-1), widths);
  1251.   return hnu0-dhnu + widths(cum);
  1252. }
  1253.  
  1254. /* ------------------------------------------------------------------------ */
  1255. /* random worker routines */
  1256.  
  1257. func multi_no_dups(x)
  1258. /* DOCUMENT xnd= multi_no_dups(x)
  1259.      returns its input vector X with any duplicate values removed.
  1260.      X must be non-decreasing and of length at least two.
  1261.  */
  1262. {
  1263.   return x(grow([1], where(x(dif)>0.0)+1));
  1264. }
  1265.  
  1266. func _multi_times(tmaster, tslave, terror=)
  1267. /* xxDOCUMENT times_list= _multi_times(tmaster, tslave)
  1268.      return the subset of TMASTER which occurs in TSLAVE.  Times in
  1269.      TSLAVE which do not appear in TMASTER are ignored.
  1270.      The optional TERROR keyword determines the accuracy with which
  1271.      TSLAVE must match TMASTER in order to be included.  TERROR may
  1272.      be a scalar, or an array of length TMASTER.  By default TERROR
  1273.      is 0.01 of the spacing between successive elements of TMASTER
  1274.      (see code for the exact expression).
  1275.  */
  1276. {
  1277.   /* handle single points and other common annoying cases */
  1278.   if (numberof(tslave)<1) return [];
  1279.   if (dimsof(tslave)(1)!=1) tslave= tslave(*);
  1280.   if (numberof(tslave)>1 && anyof(tslave(dif)<=0.))
  1281.     tslave= multi_no_dups(tslave(sort(tslave)));
  1282.   if (dimsof(tmaster)(1)!=1) tmaster= tmaster(*);
  1283.   if (numberof(tmaster)<2) {
  1284.     if (numberof(tmaster)<1) return [];
  1285.     if (is_void(terror)) terror= max(1.e-6, 1.e-6*abs(tmaster));
  1286.     if (anyof(abs(tslave-tmaster)<=terror)) return tmaster;
  1287.     else return [];
  1288.   }
  1289.  
  1290.   /* create default terror widths if not supplied */
  1291.   if (is_void(terror)) {
  1292.     dt= tmaster(pcen)(dif);
  1293.     dt(1)*= 2.0;
  1294.     dt(0)*= 2.0;  /* never same as dt(1) since at least two times */
  1295.     terror= 0.01*dt;
  1296.   }
  1297.  
  1298.   /* get list of indices of largest and smallest possible match
  1299.      of tslave in tmaster --
  1300.      the required points are where these are equal */
  1301.   list1= digitize(tslave, tmaster+terror);
  1302.   list2= digitize(tslave, tmaster-terror)-1;
  1303.   list= where(list1==list2);
  1304.   if (numberof(list)) return tmaster(list1(list));
  1305.   else return [];
  1306. }
  1307.  
  1308. /* ------------------------------------------------------------------------ */
  1309.